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ABSTRACT The structure of infrared cirrus clouds is analyzed with 
Laplacian pyramid transforms, a form of non-orthogonal wavelets. Pyra- 
mid and wavelet transforms provide a means to decompose images into 
their spatial frequency components such that all spatial scales are treated 
in an equivalent manner. The multiscale transform analysis is applied to 
IRAS 100 fin l maps of cirrus emission in the north Galactic pole region 
to extract features on different scales. In the maps we identify filaments, 
fragments and clumps by separating all connected regions. These struc- 
tures are analyzed with respect to their Hausdorfl dimension for evidence 
of the scaling relationships in the cirrus clouds. 


INTRODUCTION 

In general, interstellar clouds appear inhomogeneous with features on different 
scale size depending on the spatial resolution of the telescope and the tracer used 
to probe the structure. This structure results from a fragmentation process that 
is not well understood and cannot be easily modeled because of the non-linear 
nature of the hydrodynamic equations which describe the cloud evolution. The 
formation and evolution of interstellar clouds is governed by a variety of forces 
including gravity, magnetic fields, rotation, thermal pressure, turbulence, and 
systematic motions. In addition discrete stellar sources inject energy into the 
clouds in the form of winds and shocks, producing systematic and turbulent 
motion. A highly “clumpy” or fragmented structure may reflect gravitational 
fragmentation or the presence of turbulence resulting from the cascade and redis- 
tribution of energy injected on different scales. In contrast, a highly lilamentary 
structure might result from ordered magnetic fields. 

To study the structure and dynamics of interstellar clouds observers have 
mapped them on different scales and to different degrees, mainly with the iso- 
topes of CO, CS, and the 100 and 00 ytm infrared dust emission, each of which 
traces a different aspect of the density and mass distribution. Interpreting these 
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maps has provon to be more complicated because one does not know a priori 
what kind of features to extract and exactly how the distribution of features, 
their sizes and number, relate to the dynamics of the cloud. 

To interpret the maps one needs some method to extract structural features 
anil determine their number, size, distribution, mass, and energy content. Sev- 
eral approaches have been applied and discussed in the literature including a 
search for connected objects (Dickman, Horvath, and Margulis 1990; Falgarone, 
Phillips, and Walker 1991), the autocorrelation lunction, the structure function 
(Kleiner and Dickman 1984 and 1985), an iterative fitting procedure assuming 
Gaussian-shaped clumps (Stutzki and Gusten 1990) and structure tree analysis 
(Houlahan and Scalo 1990 and 1992). Recently we suggested that Laplacian 
pyramid (multiscale) transforms (LPT) are more suitable for determining the 
structure of interstellar clouds (Langer, Wilson, and Anderson 1993). 

The best known example of the multiscale transform is the wavelet trans- 
form (Crossman and Morlet 1987; see also the reviews by Daubechies 1992; 
Fargo 1992). Multiscale transforms provide a mathematically consistent way 
to extract structural components and map properties from astronomical im- 
ages. 1 he multiscale transform has been characterized as a generalization of 
the Fourier Transform which is capable of representing a function in terms of 
spatial and frequency localization. The localized structures at different scales 
are easier to identify in the transformed space than in the original (a;, y) space, 
llie wavelet transform is well suited to provide detailed information and deep 
insight into structure. Orthonormal wavelets have been used to extract particu- 
lar features from complex astronomical data, including: galaxy counts (Slezak, 
Bijaoui, and Mars 1990; Martinez, Paredes, and Saar 1993); stellar photometry 
in globular clusters (Auriere and Coupinot 1989); 13 CO spectral data of the 
L1551 outflow (Gill and Henricksen 1990); and, photometric analysis of galaxies 
(Coupinot et al. 1992). 

Here we apply a multiscale transform analysis of infrared cirrus clouds in 
the North Polar region using the IRAS Sky Survey Atlas (Wheelock et. al. 
1994) plates and analyze the fractal structure and morphology of the clouds. 
Infrared 00 and 100 pm IRAS maps have been analyzed for scale-dependent 
morphology by a number of authors using algorithms to characterize the feat ures 
which rely on connecting contours in the original (x, y) intensity maps. Bazell 
and Desert (1988) analyzed the fractal structure of interstellar cirrus using the 
Sky flux plates for three regions: one above and one below the plane ( (> = -f23° 
and —14°) and a third region at b = —40° containing two well known high- 
latitude MBM clouds. Dickman, Horvath, and Margulis (1990) analyzed IRAS 
images of five molecular cloud complexes (including the well known regions of 
Chameleon, p Oph, and Taurus). Both studies found a highly fractal structure 
for the clouds and concluded that the clouds had a highly turbulent structure. 
However, they differ in the behavior of the scaling from region to region, which 
might not be surprising considering that one concentrates on molecular and the 
other diffuse regions. Our multiscale analysis of cirrus emission reveals a more 
complex structure than either of these previous studies. 
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MULTISCALE TRANSFORMS 

First, we review the various multiscale transforms (Pyramids and Wavelets) 
and compare them to Fourier Transforms. Laplacian pyramid transforms pre- 
ceded and spurred the recent interest by mathematicians in the more formalized 
orthogonal wavelets. In practice the orthogonal wavelets have proven useful pri- 
marily in image data compression, while the Laplacian pyramid has remained 
more useful for carrying out a variety of image analysis tasks. 

Fourier Transform 

A Fourier Transform decomposes a function f{x) into a linear composition of 
Fourier vectors, whose basis functions are sines and cosines, defined by their 
Fourier coefficients, 


f(k) = (2n)- l/2 f f{x)e ik *dx. (1) 

Unfortunately t ikx oscillates forever and the information content of f(x) is de- 
localized among all the spectral coefficients f(k). Regions of f(x ) that vary 
sharply - for example a delta function - have their information spread among 
al l of k space. Such transform properties are not very useful for real astronomical 
images which contain a great deal of structure. 

One approach to achieve localization is to use a window to isolate the por- 
tion of the frequency or spatial scale of interest. Windowed transforms can also 
be applied sequentially to recover all of the function, however such transforms 
use fixed window size and thus do not analyze the function at all scales. An 
example of a windowed transform is the Gabor transform, 

G a (x) = (4ira)-V 2 e (2) 

which uses a Gaussian window. The real part ol the transform of the Gabor 
function, lieG w is shown in Figure 1 for three values of u>; note that the Gabor 
transform has a fixed window size and that the transformed G oscillate within 
a fixed envelope. 

Wavelet Transform 

The wavelet transform localizes information in space and scale and the decom- 
position is done by translation and dilation of a single “Parent ’ function, the 
analyzing wavelet. The wavelet transform is defined by 

T tJ {a,b) — |a| -1 / 2 / f(x)g{(x - b)/a}dx (3) 

where the factor a covers different ranges in scale (frequency) and b moves the 
localization center (position or time). The wavelet transform separates data or 
functions into different spatial components and analyzes each component with 
a resolution to match its scale. The wavelet transform function must have the 
property that its mean is zero over .t, 

/ (j(x)dx — 0. 


( 4 ) 



R«G' ob (x) 
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FIGURE 1 The functions # a 6 and RtG^t for three values of a and u> ( b = 
0). The Gabor transform has a fixed window size in contrast to the wavelet 
transform’s flexible window (Martinez et al. 1993). 


One of the best known examples of a wavelet transform in astronomy is the Marr 
(or “Mexican Hat”) wavelet, which in one dimension has the form, 

g(z) = (1 - z 2 )e~ z ' 2 /' 2 , where z = (x - b)ja (5) 

The transformed function <j a t> is shown in Figure 1 for three values of a (assuming 
6 = 0) where one can see that the Marr wavelet has a flexible window size. 
Wavelets have the useful property that they preserve scaling behavior and are 
sensitive to signal variations but not to constant behavior. They also provide 
a useful description of turbulence because they retain information about the 
spatial structure of the How (c.f. the review by Farge 1992). 

Discrete Wavelet Transform 

I here is a discrete wavelet transform (DWT) analogous to the discrete Fourier 
transform (DFT), both of which are rotations in function space, but take an 
input space into different domains. For the DFT the basis functions are unit 
vectors or Dirac delta functions (in the continuum limit). For the DFT the 
domain has sine and cosine basis functions, while for the DWT they are localized 
in space. The best known examples of the DWT use the Daubechies basis 
functions (Daubechies 1992). We do not have space to review the DWT (see 
Press et al. 1993) but point out that the transform is a matrix whose odd and 
even rows perform different convolutions: the odd rows generate components of 
the data convolved with filter coefficients and perform a moving average with 
the smoothing filter H; the even rows perform a convolution with a filter L 
which represents the detailed information in the data. The DWT is applied 
hierarchically to the data vector (i.e., first to the full length data N, then to 
the “smooth’ vector ol length N/2, etc.) ultimately producing a set of wavelet 
coefficients that represent the smoothed information on the largest scale and 
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FIGURE 2 Scale invariant filter bank: Filter responses, F n (k) — Fq( 2 n &), on 
the left. 

detailed information on scales differing by factors of two. Specific examples are 
given by Daubechies (1992) and Press ct al. (1993). 

Pyramid Transforms 

Pyramid transforms preceded and spurred the recent interest by mathematicians 
in the more formalized orthogonal wavelets. However, one of the major difficul- 
ties with the orthogonal wavelets is that, while their basis functions display shift 
and scale invariant properties, the coefficients in these expansions do not (see 
Strang 1989; Simoncelli et al. 1992). The fact that the power within a given 
scale is not invariant to translations of the input should be enough to make one 
wary. In addition, when implemented as separable filters, the orthogonal wavelet 
transform creates three highly anisotropic bandpass components at each scale, 
which have no simple relationship to underlying physical properties. 

Pyramid transforms provide a means to decompose images into their spatial 
frequency components such that all spatial scales are treated in an equivalent 
manner. The concepts are equally applicable to data of any dimension, where 
each axis represents some continuous physical parameter such as spectral wave- 
length, time, or in the case of images vertical and horizontal lengths or angular 
displacements. An intuitive grasp of the idea can be obtained in one dimension 
by considering the illustration in Figure 2. The horizontal axis represents a pa- 
rameter such as space and the vertical axis spatial frequency. If we consider a 
signal 5(tf) that has been bandpass limited by a fixed amount, then it can rep- 
resented by a uniform sequence of samples S(xi), as set by the Nyquist sampling 
theorem, with no loss of information. Now consider passing the signal through a 
bank of filters that are scaled copies of one another, such that F n (k) = Fo{2 n k), 
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where the scaling by factors of two is chosen for convenience. The output of 
each filter treats the signal in a scale invariant fashion because of the scaling 
relationship between the filters. The outputs of each filter must be sampled at 
a rate that is proportional to the corresponding bandwidth, hence each band is 
sampled at a rate that differs from its neighbors by a factor of two, as illustrated 
in Figure 2. This method of signal decomposition is the basis of all scale invari- 
ant multiresolution representations such as Pyramids and Wavelets. The signal 
processing community have labeled them as a special class of subband encoding 
schemes. 

The generality of this form of signal analysis is best illustrated by noting 
that the human visual and auditory systems utilize it. The scale invariant nature 
of many physical processes is what makes these representations important. The 
statistics of the signals and the forms of the structures produced by these filters 
provide measures of how underlying physical processes change with scale. The 
coefficients in these representations are much more statistically independent of 
one another than they are in the original data formats, which leads to a rich set 
of localized descriptors of images and better data compression for storage. 

The number of such decompositions is large since there are many filter 
designs whose scaled filter banks will cover the frequency range of interest with 
sufficient density to prevent loss of information. This means there is a multitude 
of possible multiscale transforms unlike the Fourier transform. The choice of 
which filter to use is determined by factors such as efficiency of computation, 
signal-to-noise ratios, information storage and the type of data analysis one 
desires to perform. When one goes into higher dimensional spaces the choices for 
the shape of the filter become increasingly more flexible and hence increase the 
number of possible transforms. The orthonormal wavelets, which have recently 
received a lot of recognition, are a particular subset of the multiscale transforms. 
These provide critically sampled representations with minimal storage, as well as 
having some nice mathematical properties, but the constraint of orthogonality 
leads to filter designs that are not necessarily the best for many applications. 
There are a number of overcomplete transforms, such as the Laplacian pyramid, 
based on circularly symmetric filters, oriented pyramids, which have “wavelet 11 - 
like filters, and a recent new class with a property called “shiftability” (cf. Strang 
1989, Simoncelli et al. 1992), 

The Laplacian pyramids are best understood as the outputs of a bank of 
scaled circularly symmetric, bandpass spatial frequency filters, with center fre- 
quencies k n — A,'o/2 n and roughly octave bandwidths. The simple variant used 
here is called an FSD (Filter, Subtract, and Decimate) pyramid (cf. Van der Wal 
1991). Starting with the original image, designated as G'o(x, v/), the following 
rules are applied recursively to create a sequence of lowpass images (or Gaussian 
levels) Ci H (x, y) and bandpass (or Laplacian levels) L n (x,y): 


G w + 1 — H * G n {Filter) 

(6) 

L n — G n - G ' n + 1 {Subtract} 

(7) 

G n + 1 = {Decimate} G n + 1 

(8) 
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The filter operation H * G n involves convolving the image G n with a lowpass 
filter H . A separable filter, H(x J y) = h(x)h(y), was used in this work, where 
we adopted a five tap filter for h(x) having the tap values 1/16, 1/4, 3/8, 1/4, 
1/16, which produces an approximately circularly symmetric filter in the spatial 
frequency domain. The Laplacian components, X n , are computed by subtract- 
ing the low pass version from the unblurred one at each scale. This operation 
is equivalent to filtering with a “Mexican hat 11 or difference of a Gaussian-like 
shaped kernal. The blurred version, G n , is then subsampled by throwing away 
every other pixel and row, which is traditionally called “decimation.” Decima- 
tion is justified because the lowpass filter reduces the spatial frequency content 
such that little aliasing is introduced by this process. Since this reduces the 
number of pixel values by four at each stage in a two-dimensional map, the com- 
putation of all the levels L\ to Ln takes only one third the time and resources 
to compute Lq. Typically the final level N is set by stopping the process when 
the smallest dimension of the array G n +\(x, y) would be no smaller than eight. 
Mathematically, the set of values L n (x,y) for n = 0 to N plus constitute 

the components of an overcomplete non-orthogonal wavelet decomposition of 
the original image G 0 . The bandpass components show less than two percent 
anisotropic response to spatial frequency orientation and less than five percent 
to translation of the input image. This operation is a computationally efficient , 
basic wavelet transform that has a relationship to some simple multiscale statis- 
tics. The G n are mean values computed over regions that increase in size by 
factors of two and the L n provide the local deviations from these multiscale 
mean values. 


IRAS IMAGES OF THE NORTH POLAR CIRRUS 

To analyze the structure of cirrus emission we have chosen a region in the north 
Galactic pole (declination 90° at map center) in order to minimize the line- 
of-sight confusion. We obtained 60 and 100 /an images of the North Polar 
Cirrus from the ISSA data measuring 5° on a side with pixels of 1.5 / xl.5 / , 
although the actual resolution of the images is closer to 4' (the IRAS beam size at 
100 /mi). Figure 3 is a contour plot of the 100 /an image and shows the highly 
complex filamentary and globular nature of the emission. From the 100 and 60 
/mi images one can calculate a dust color temperature and opacity assuming 
a one component model of the dust grains (cf. Langer et al. 1989). The 
color temperature map is roughly uniform with a temperature of 25 K arid the 
temperature drops in the regions of brightest emission where the opacity is a 
minimum. This behavior is well known for the 60 and 100 /mi maps and is a 
result of using an overly simplified grain emission model. Therefore we decided 
to restrict our analysis to the 100 /mi maps as the 100 //m emission arises from 
cooler dust, is optically thin, and thus more likely to arise uniformly from all the 
material along the line-of-sight. The 100 /mi maps show a great deal of structure 
on all scales and the variation in emission is either due to increased column 
density, changing dust grain distribution, or variations in heating sources. 
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FIGURE 3 The 100 /an image of the North Polar Cirrus (dec = 90°, 0^ is 
at the top and RA goes counterclockwise with 6 ^ on the left); the axes are in 
arcmin with respect to the center. Intensity units range from 2 to 16 x lO 6 . 

RESULTS 

The LPT applied to the 100 /mi image of the North Polar Cirrus generates a 
set of amplitude maps for each of the basis functions, ranging from L 0 to L 5 
covering the smallest to largest scales (1,5' to 9'), and G 5 (the original image 
is labeled G 0 ). The G n are mean values computed over regions that increase 
in size by factors of two and the L n provide the local deviations from these 
multiscale values. The L 11 images represent the detailed information and have 
equal positive and negative areas (average over the maps is zero). 65 , which has 
a resolution of 48', is smooth and positive definite over the map. There is not 
room here to display the Laplacian images but examples of the transform maps 
for CO clouds can be seen in Langer et al. (1993). 

To extract the structure in the map we assume that at each scale the features 
are embedded one within the other, or superimposed along the line-of- sight. We 
define a clump or filament as a connected region of positive amplitude, that is 
the zero contour boundary in the Laplacian maps. Bubbles or valleys would 
correspond to connected negative regions. (These maps are better for isolating 
features because of the separation of different scales and the preservation of 
locality in the transformed space.) 

For each image we can determine the following properties of the connected 
positive or negative regions: 
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Perimeter = I\ 

Area — A, 

Intensity = f</>(x,y)dA, 

where </>(x,y) — amplitude of the wavelet basis function at (x,y), 

Brightness = [J </>(x, y)d A]/ A, 
rms - J\</>{x,y)\ 2 dA. 

The major problem with this approach is that there are many features at 
small scale that correspond to noise in the maps or are on smaller scales than the 
beam resolution (the pixel area is 2/25 arcmin 2 while the IRAS beam is roughly 
3' to 4'). We avoid this problem by only considering features whose area, /l, is 
greater than l(i arcmin 2 ). Table 1 summarizes the number of distinct features 
of positive amplitude for each of the Laplacian maps. 


TABLE 1 Summary of Cirrus Features in Each LPT Basis Map 


Order 

(i) 

(AA)* 

(') 

Number of 
Features 

0 

1.5 

1423 

1 

3.0 

084 

2 

0.0 

252 

3 

12.0 

81 

4 

24.0 

28 

5 

48.0 

11 


The wavelength of each band is approximated by the width. AA, ol each high 
pass filter. 

Here we concentrate on only one property of the structure analysis of the 
cirrus clouds, the scale-dependent morphology, or fractal structure. One widely 
used approach to measuring the morphology of the cloud structures is the Haus- 
dorfr dimension (Bazell and Desert 1988, Dickman et al. 1990) which is the 
fractal dimension of the objects in the cloud. The fractals are self-similar under 
scale changes and so a determination of the Hausdorff dimension, D, can help 
identify the scale-dependent morphology of the cloud. Ideally, the geometrical 
structures and the scale dependence provide information about the forces at 
work in the cloud. D is defined as 

A 1 /' 2 = K1 a/d (9) 

where a plot of log perimeter versus log area yields a slope ol D /2 and has 
an intercept equal to -D(loyK). Regular geometrical objects (circles, ellipses, 
squares, etc.) all have D - 1 (i.e., the same scaling) but different intercepts 
(i.e., K depends on shape). A value of D « 4/3 characterizes the relationship 
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expected for Kolmogorov turbulence (incompressible, homogeneous, isotropic 
turbulence); see the discussion by Dickman et al. Filamentary structures that 
scale only with length, on the other hand, have I) & 2, as would be characteristic 
of gas supported by magnetic fields in the radial direction. Figure 4 shows an 
example of the plot of log(P) versus log(A) for the / = 1 image ( AA « 3'). Note 
that the fit is not particularly good for large areas (A > 400 arcmin 2 ). For 
the small wavelengths (/ = 0 and 1) the fits are dominated by the many small 
features in the map (note each cross in Figure 4 can represent many clumps of 
the same area and perimeter). Figure 5a plots D as a function of the l value and 
shows relatively little variation with wavelength. The log(P)-log(A) fits for / = 4 
and 5 do not have much area dependence and the average value for D in Figure 5 
is about 1.45. Bazell and Desert found a rather different average value, D « 1.25, 
for the cirrus clouds, however, there was variation within each image (see their 
Figure 2). Dickman et al.’s study of 100 /xm maps of molecular regions found 
D ranging from 1.2 to 1.3, similar to those found for cirrus clouds (remember 
that these results use a completely different method than that discussed here to 
determine the map features). 

Determining the Hausdorff dimension as a function of area, in the sense of 
considering several ranges of area within an l map, yields somewhat different 
results. Figure 5b plots D for the / = 1 image divided up into four different 
ranges. For the smallest features the Hausdorff dimension is 1.25, similar to the 
result of Bazell and Desert, but it increases with area reaching a value about 2 
for the largest objects in the image. As the / — 1 map corresponds to features 
that respond to a multiscale filter of about 3' (in at least one direction) the 
largest features resemble thin filamentary or web-like structures. This result 
suggests that the forces that control the structure at large and small scales in 
the cirrus map of the north pole region are different. 

DISCUSSION 

D for all the North Polar Cirrus features ranges from 1.35 to 1.5 with an average 
of 1.45 for all the Laplacian maps. These values are larger than the average 
value of 1.25 found by Bazell and Desert in three cirrus maps, their maximum 
value of 1.40 for plate 2, and the values for D(& 1.18 - 1.28) for 100 /xm IRAS 
maps of the molecular clouds studied by Dickman et al. It is not clear if the 
differences are due to the choice of objects or the method to extract the map 
features. We will analyze the Chameleon region in the near future to test the 
latter possibility. The fractal nature of the cloud B5 studied with CO maps 
(Langer et al.) has D ranging from 1.3 to 1.7. These differences may occur 
because CO emission probes the dense interior of the cloud where the dynamics 
are much more influenced by gravity and embedded sources. D for cirrus clouds 
is consistent with a turbulent gas (Kolmogorov turbulence including the effects 
of intermit tency — see Falgarone et al.). However, the strong dependence of D 
on the area range (Figure 5b) shows that for / = 1 and 0 (the smallest scales) 
D increases with increasing area reaching values about 2. This dependence 
suggests that the smallest features obey a Kolmogorov scaling law (D is close 
to the Bazell and Desert value) while the largest feature at / = 0 and 1 are 
filamentary-like structures. Visual inspection shows these features to consist 
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0.5 1.0 1.5 2.0 2.5 3.0 3.5 


Log (A) 

FIGURE 4 Determination of the Hausdorff dimension, D , and the shape fac- 
tor, A', from plots of perimeter versus area for the / — 1 image. The fits are 
listed in the Figure. 


2.0 
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1.2 
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Wavelength 

FIGURE 5a. Hausdorff dimension, D, FIGURE 5b. HausdorlT dimension, D. 
versus wavelength, see Table 1) versus different ranges in area for the 

using all the features in each map. I — 1 map. 
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of long single filaments and web-like structures. The latter are likely due to 
the overlap of separate features along the hue of sight since the data have no 
velocity information and/or the merging of very many small scale features within 
the filters. Higher resolution 100 pm data is needed to test this hypothesis. 
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